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ABSTRACT 


This thesis is a partial analysis of the Naval Space Command statistical 
orbit determination algorithms. Through a process called Differential Correction, 
data from space surveillance radar observation stations is synthesized with 
previously accumulated element sets to maintain accurate orbital object position 
information. Differential Correction is a nonlinear least squares process employing 
statistical techniques to minimize the residual measurement error thereby 
increasing relative position information accuracy. This study focuses specifically 
on the algorithmic methods of solution to the systems of normal equations 
generated by the Differential Correction process. A comparison and analysis of the 
present Naval Space Command method and the singular value decomposition 
method is presented. Algorithmic constructions are presented for both methods 
and problematic areas are highlighted. The principal focus herein is to demonstrate 
the benefit of singular value decomposition when attempting to solve systems of 
equations whose coefficient matrices are dense and nearly singular. These results 
generalize to commonly employed normal equation solution algorithms and are 
intended for further study and possible incorporation by Naval Space Command as 


part of future modernization plans 
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I. INTRODUCTION 


The Naval Space Command maintains a database of element sets for more than 


- 9,000 objects. Through a process called Differential Correction, data from space 


surveillance observation stations is synthesized with previously accumulated element sets. 
Differential Correction is a nonlinear least squares process employing statistical techniques 
which provide for accurate orbit state estimation. Radar measurements of an object’s 
motion are collected by dispersed observation stations and passed to NAVSPACECOM 
for central processing. NAVSPACECOM receives daily approximately 270,000 new 
observation sets and uses them to update as many as 18,000 element sets. 
NAVSPACECOM performance reports indicate that approximately 98.5% of the element 
sets get updated without manual intervention by their computer software called 
AUTODC. The remaining 1.5% must be manually analyzed. Differential correction is a 
highly complex step-wise process that entails more than 16 separate mathematical 


computations. Even with automated support, this constitutes a significant work load. 


Numerical linear algebra is at the core of the differential correction process. In 
general, data is accumulated, normal equations are formed, and numerical linear algebra 
routines are called to solve the matrix system of normal equations. The method of normal 
equation solution and its associated algorithms are the focus of this analysis. The purpose 
of this thesis is to demonstrate the benefits of using singular value decomposition when 


obtaining least squares solution to systems of normal equations. 














Ill. BACKGROUND 


A. STATISTICAL ORBIT DETERMINATION 


The NAVSPACECOM differential correction method is a sequential batch 
nonlinear least squares statistical process. This is an iterative method which refines the 
stored orbital element sets by applying state adjustments obtained from differential 
correction of the current orbital observations. These state updates are tested for fit in the 
least squares sense and if acceptable applied to the previous nominal state then stored as 
the new nominal orbital element set. The fundamental mathematical steps are now 


outlined; for further details, see Vallado [1] and Danielson and Canright [2]. 
B. DERIVATION OF THE NORMAL EQUATIONS 


These descriptions have been simplified and are only intended to provide sufficient 
background for the purpose of detailing the algorithmic composition and term-wise 
structure of the normal equation system. For a very detailed description of the entire set 


of NAVSPACECOM procedures refer to Danielson and Canright [2]. 
C. GENERAL LINEAR LEAST SQUARES 


We begin the process by noting that our goal is to solve the generally inconsistent 


system of equations whose matrix representation is 


Ax =b (1) 











by minimizing the residual, r, where 
r=b-—Ax (2) 


The method of least squares can be applied to solve such linear systems. Typically 
these systems are highly overdetermined. They have many more rows (m), or equations 
than columns (n) or unknowns and as such are inconsistent. From elementary linear 


algebra [3], we know that the range of the matrix R(A) is the orthogonal complement of 


the null space of its transpose N( A’ ). As such, the multiplication A’ always produces a 
consistent set of n equations and n unknowns. The only issue with this method however 
arises when some set of the variables from the original linear system are correlated. If this 


happens, the newly formed consistent system will fail to have a unique solution. 
The general method proceeds as follows. Left multiplication of equation (1) yields, 

A™Ax = A’b (3) 
Equation (2) implies 

A’r = A’(b-— Ax) = 0 (4) 


Equations from the original system may be weighted by applying a diagonal weighting 


matrix. This expands as 
A™WAx = A‘*Wb (5) 


Residuals from the new, weighted system satisfy 
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A’Wr = A'W(b- Ax) = 0 —. (6) 
Provided A’ WA is nonsingular, we obtain the solution as, 
x =(A7WA)" A7Wb (7) 


This result commonly referred to as the normal equations forms the basis of the least 


squares solution process. 


Differential Correction employs sequential batch techniques when acquiring the 
current state solution from different observational data input sets. This method 


synthesizes the separate batch solutions as follows. 
Given the systems for two batches as 
A,.x = b, and A,x = b, 


In order for the single solution to be meaningful, these batch systems must be 
simultaneously solved. This is accomplished by forming the weighted least squares system 


as follows. 
[(A7 WA), +(A7WA), x = (A? Wb), +(A’ Wb), 
The process can be simplified slightly by reusing the solution from x, 


[(A7 WA), +(A7WA), x = (A’WA),x, +(A’ Wb), 








When acquiring the weighted solution simultaneously, different weights may be applied to 
the different sets of equations either by reducing W, or by scaling (A’ WA), . Notice here 


that the batching process does not significantly alter the basic least squares problem nor 


does it alter the complication arising when some portion of the left hand side constitutes a 


singular matrix. 
D. LEAST SQUARES APPLIED TO OBITAL MECHANICS 


The reader is directed to Danielson and Canright [2] for extensive mathematical 
description and existing code documentation from which the following outline extract was 


taken. 


(i) Assume an initial nominal state 
nominal ~ 


(ii) Compute the values of the observed parameters Y, at N times corresponding to the 


observations Y, (each observation set contains at a maximum 6 numbers) 


Q7 c) 


Y= 43 , ¥.=] : | , where7=1,...N 


z 


(iii) Compute the residuals or “O-Cs” (Y,-—Y.); (observed minus calculated 


parameters). Arrange these as the 6Nx1 column matrix 








(iv) Calculate the partial derivatives 


eY, __aY, Alr,v) 
aX a(r,v) ax 























On oY; 

Oy; 

oY, OY, OY, OY, OY, OY, |; % Os 
Or; Or, Oy OY, ov, Vy Og | Ox 
aX, aX, 

= ware 
Oe Dee Teg Meg Me Mee OM, 
Or, Or, = Oe = OV, SOV, Oe || Ov, OY, 
“a 

Ve We 

gee 


at each observation time and arrange the coordinates of the position and velocity vectors, 


(7;,7),7x) and (v;,v,,Vx), into the following 6Nx8 matrix. 


aY,, 
aX, ); 
aY, aY, 
Ae eee k ass 2 
ax, } aX, } 








(v) Form the normal equations: | 


A*WAx = A'Wb 








Where W denotes a 6Nx6N weighting matrix. 


Finally 


are the differential corrections. Note that A’ WA is an 8x8 matrix, and that A’ Wb is an 


8x1 matrix. 


(vi) Solve the normal equations: 
x =(A’ WA)" A’ Wb 
(vii) Update the elements: 
Anew = Sag + % 
(For the first iteration, Xo is Xnominal-) 


(viii) Lastly, apply the following RMS test to determine if iterations should continue. 


Rus, - RMS. < 
RMS, | 


3 








| 1m? b b 
where, RMS; = wa = N-1 
i=] 


Generally, in nonlinear least squares our goal is to mimimize the sum of the 


residuals squared. It is therefore appropriate to use some measure of this as stopping 
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criteria. In the differential correction process, simply cease iteration when the RMS 
“stops” changing. Since all of the necessary input to the RMS equations has been 


accumulated this proves an efficient method as well. 
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Ill. NORMAL EQUATION SOLUTION METHODS 


A. GENERAL 


At present the Naval Space Command, Differential Correction process relies on 
Gauss—Jordan elimination with full pivoting for solution to the normal equations. The 
subroutines responsible for processing the complete normal equation solutions are named 
AMA06, AMA04 and AMAO3. After the appropriate Differential Correction algorithms 
have formed the standard system of normal equations as described 1n the previous section, 
AMAO6 is called to read in the initial array entries, form the matrix E =A’ WA, duplicate 
it and begin preconditioning processes designed to identify singular matrices. Following 
sufficient preconditioning, AMAO6 then calls AMA04 which in turn computes E™. 
Finally, AMAO3 is called to apply Gauss-Jordan elimination and solve, 


x =(A’WA)'A’ Wb. The 8x8 matrix A’ WA is input as E and the 8x] matrix A’ Wb 


is input as G.” 


The following section describes specific processes and algorithms used to solve the 
normal equations with Gaussian elimination and with Singular Value Decomposition, 
SVD. Each process is initiated with the previous normal equation derivation. The present 
NAVSPACECOM differential correction process uses the Gauss Jordan elimination 


method with full pivoting. 
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Focusing back on step (vi) of the least squares orbital mechanics process, solution 
of the normal equations, we now compare the present NAVSPACECOM method with 
Singular Value Decomposition method. The refined problem is now, how in acquire the 
best approximation at each step of the nonlinear approximation process so as to ensure 
maximum efficiency of the iterative routine while minimizing the associated rounding 


error in the next iterate? 
B. NAVSPACECOM GAUSS-JORDAN ELIMINATION 


The Naval Space Command uses Gauss-Jordan elimination with full pivoting for 
solution to the system of normal equations. Their algorithm is very similar to the one in 


the book, Numerical Recipes, by Press, et al. [7]. 
Typical Gauss-Jordan with Full Pivoting 


Every Gauss-Jordan step is a left multiplication by an elementary matrix, [3]. An 


overview of the general method of Gauss-Jordan elimination is as follows: 
(1) Augment the right hand side with the nxn identity matrix. 


(ii) | Perform elementary operations on the augmented, nx2n system until the 


nxn identity matrix is reformed in the first n columns of the matrix 


(iii) Recall the last n columns of the augmented matrix as the inverse. 
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This method 1s usually done in place; that is, the actual augmentation is never 
performed. The algorithm simply replaces the original input matrix, in place, with its 
inverse. For this reason, if there is ever a need to recall the original data, a copy of the 


original matrix must be made on input. 


Full pivoting is the process by which rows and columns of the original matrix are 
reordered so that the element largest in magnitude is moved to the or left corner of the 
matrix. Then, each successive pivot row in maneuvered similarly. Why do this? In the 
computational process, the subsequent rows below the pivot in question are being reduced 
to zero by dividing through by a scaled multiple of that pivot. If the process is not begun 
in this manner, that is to say if the element below the pivot is greater in magnitude then the 
pivot, the pivot must be scaled by increasing its magnitude. The scaling factor is the 
element being zeroed out divided by the pivot, at if the pivot is small enough this 
| approaches division by zero. In finite precision arithmetic, this can lead to a divide by zero 
condition. If the difierencs in magnitude of the pivot is sufficiently less than the element 
being zeroed out, the computer can evaluate the expression and ee “Not-A-Number”. 
Even if the division by zero extreme condition does not occur, scaling the matrix through 
by these sufficiently small quantities inflates the rounding error. How does pivoting 
protect against this? By always scaling so that each successive pivot is the next largest in 
magnitude, the previous process of dividing the subdiagonal elements, results in scaling 
the pivot row by decreasing its magnitude. This has the effect of moving the division 
operations away form the 1/zero condition and deflates the magnitude of rounding error to 
the minimum possible given the actual element-wise composition of the matrix. 
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Partial pivoting is a similar process that allows only row exchanges. While this 
may appear to be a sound concept with less computational overhead, in practice, matrix 
elements, not directly below (in some other column) the pivot being operated a may still 
be sufficiently greater in value than subsequently available pivot choices. If this occurs, 
the solution is again driven toward a divide by zero condition. By manipulating both rows 
and columns, as in full pivoting, we reduce the likelihood of the divide by zero condition 


to the greatest extent possible. 


When employing the full pivoting routines, the general Gauss-Jordan process is 
augmented by left and right permutation matrices which track the pivot maneuvers. The 
corresponding information is stored as an array and applied during the back solve process 
thereby ensuring correct association of the coefficients and variables from the original 


system of equations. 


It must be noted that the previous pivoting processes only reduce the likelihood of 
algorithm failure. They do not prevent it. In the event of near singular, noninvertible 
conditions, the Gauss-Jordan full pivoting routine has no mechanism to correct or 
compensate. This is an inherent critical deficiency. Failure due to nearly singular input 
matrices routinely occurs when processing least squares systems that are derived from 


closely sampled data. 
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Unique features of NAVSPACECOM Gauss-Jordan 


The following system processes are extracted as reported in the code 
documentation research of Danielson and Canright, [2]. The matrix representation of the 
system of normal equations is loaded as the input matrix FE. The input matrix is symmetric 
positive semi-definite by its inherent construction. As such, only values on and to the right 
of the diagonal are read as input. The below diagonal elements are copied from their 
respective symmetric counterparts. At this point a preconditioning routine, AMAO6, is 


called. It samples the matrix to determine if any off-diagonal element of E is too large 


relative to its corresponding diagonal elements. AMAO6 tests each E; > (SING)E:E;, 
where SING is a parameter value. If any E; violates this inequality then that row (for 


even-number calls) or column (for odd-number calls) is “inactivated”. The threshold 
SING 1s initially set to near one. Then, if the active matrix is singular, the threshold SING 
is lowered by 10% (attempting to inactivate more rows/columns) and the solution is tried 
again. At this point, a saved copy of the original matrix must be reloaded. The immediate 
effect of this process is to eliminate divide by zero conditions when pivot elements are 
very small. Eventually, either a solution is found, or the whole matrix is made inactive 
(flagged by a return value SING=0), or the threshold gets too small (SING < 0.01) 
indicating that the entire system is singular. In the event that no matrix remains, an error 
condition is returned and all further processing of this set of observational data ceases. A 


tally of these failures is presented as an output report in the batch run statistics. 
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There is a fundamental problem with this preconditioning methodology. Under no 
circumstances are the root causes of the perceived singularities dealt with. Further, in 
terms of computational effort, the computer is permitted to cycle excessively without 
sufficient intermediate checks to determine if the initial perceived singular conditions are 


removable. 


In summary, Gaussian elimination is generally the most computationally efficient 
numerical method available; however, it suffers: from numerical instability under certain 
conditions, which are detailed in Chapter IV, Comparison of Methods. Unfortunately, 
these very same conditions generally arise in systems of normal equations. As noted 
earlier, whenever correlation of variables from the original system equations occurs, the 
resulting matrix A7A represents a consistent set of equations with redundant solutions. 


These redundant solutions drive the matrix toward a singular condition. 


C. SINGULAR VALUE DECOMPOSITION 

The Singular Value Decomposition, SVD of an arbitrary mxn matrix is the 
factorization of A into UXV", where U and V are orthogonal matrices and 2 is the 
rectangular mxn matrix whose first r rows form a square, diagonal submatrix with 
elements o,---o, , i.e., the singular elements of A with the remainder of & being zeros. 


The process proceeds generally as follows (See [6]): 


(1) Reduce the general input matrix A to bidiagonal form B with orthogonal 


matrices U and V where A = UBV’ . B is nonzero only on its main 
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diagonal and first super diagonal. This is accomplished with Householder 


transformations in the actual algorithms. 
(ii) Find the SVD of B, B= U SV", where © is the diagonal matrix of singular 
values and U, and Vi are orthogonal matrices whose columns are the 


respective left and nght singular vectors. The Gram-Schmidt process . 


produces this result. 
(iii) Combine these decompositions to form A= (UU uv. The 
columns of U= (U U_)and V= (V V_ )are the respective left and right 


singular vectors of A. 


Step (i) reduces the A matrix to bidiagonal form by applying Householder 





transformations on both left and right sides. The symmetry of the original matrix is 
preserved throughout. Step (11) employs the Gram-Schmidt process to orthogonalize the 


columns of U and V. Step (iii) reforms the matrix factors. Specific algorithms for 


these subordinate functions are detailed in Golub [5] and Demmel [6]. The SVD 
algorithm simply incorporates them. Another extremely useful secondary benefit of the 
SVD is provided by its fundamental structure. That is the factor matrices U and V have a 


very special structure. They are orthogonal. Where, A = ULV’ , the columns of U are 





the eigenvectors of AA” and the columns of V are the eigenvectors of A’A , Strang, 


[3]. Further, orthogonal matrices have the nicety that U7 U = UU? =I. This infers 
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directly that when solving the system of normal equations, 


Ax =b as ULV‘'x =D therespective transposes, U* and V, are left 


multiplied yielding x = V=‘U"b. Interestingly enough, no actual inverse matrices have to 
be calculated. Further, the transpose of U and V can be done in place. The real benefit is 
the complete absence of rounding error when taking the transpose of a matrix. 

The problem with the Gauss-Jordan method is the requirement for the columns of 
A to be independent. As described in detail in the previous section, when ill-conditioned 
matrices are input, a great deal of computational effort is required in preconditioning the 
matrix. Without this preconditioning, even nonsingular but very ill-conditioned matrices 
may cause algorithms to break down. The computer simply perceives the matrix to be 
singular to its level of machine precision. | | 

The key to SVD’s computational stability lies in its orthogonality. Matrix rank 
problems arise frequently in computer arithmetic. Determining the rank of an ill- 
conditioned matrix can be challenging in the presence of roundoff error and noisy data. 
The SVD allows for practical dealing with neta rank deficiency 

The following theorem is taken from Golub, [5]. See [5] for a detailed proof. The 
theorem provides the detail necessary for determining how “close” the given A matrix is to 


one of lower rank. The 2-norm here is the matrix derivation of the vector Euclidean norm 


[3] defined as follows. 


least upper bound (A)=max (at where |Al < least upper bound (A) |x| 


hal 


18 








Theorem 


Let the SVD of AER™" be given. Ifk <r = rank(A), and 
ee 
A, = > ou,y; 
i=l 
then, 
min |A-Bl = la -A) cree 
rank(B) =k 

This striking result offers a method of computation using stored values which 
specifically reveal the magnitude of the degenerate numerically singular condition of the A 
matrix. The details of this theorem indicate that the smallest singular value of A is the 2- 
norm distance of A to the set of all rank deficient matrices. Iteration is no longer required 
to determine the degree of singularity. 

Now, let’s look at an example of how to make the most out of ill-conditioned least 
Squares systems with the SVD. See Strang [3] for more on this. In general, the least 
squares problem has one very stringent requirement, the columns of the A matrix must be 
independent or the rank of A must be equal to n, the number of columns. This is often 
referred to as full rank. If not, A is not invertible then Ax = b can not determine x. As 
described earlier, any vector from the null space of A can be added to x. Now let’s 
examine what happens. There are two possible situations, either the rows of A may be 
dependent or the columns of A may be dependent. The first situation implies the system 


of equations may have no solution and the second situation implies that any solution is not 


unique. The dependent column case makes this a particularly difficult yet interesting 
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problem. As discussed earlier, when we have dependent rows the solution we seek may 
be outside the column space of A. Our course of action now becomes simply project b 
onto the column space of A. Now the greater challenge. After making that projection we 
find A has dependent columns and the solution is not unique. At this point, we must now 
employ the criteria for selecting the optimal solution and choose the one with muni 
| length. | 


Consider the following example: 


o, 0 0 0 
A = 0 o, 0 0 
0 0 0 0 
0 0 0 0 


where A is diagonal with dependent rows and columns. 

Here we see the columns all end in 0. As per case (i), the closest vector to 
b = (b,,b,,5,,b, ) is p=(8,,5,,0,0)the projection onto the column space of A. The 
magnitude of error here is b = (0,0,5,,5,)the perpendicular to the columns. The best 
solution now is attained when we solve the first two equations. Since the last two 


equations indicate 0= b, and 0=5,, the error in those equations cannot be reduced but 


the error in the first two will be zero. 


07 0 0 0 Xi by 
— 0 o, 0 0 
Ax =p is Xa | | b2 (8) 
0 0 0 Ollx,} | 0 
0 0 0 Oflx,} [0 
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Now the challenge, the dependent columns imply x is not unique! The first two 


components are © and ae but x, andx, are completely arbitrary. Now apply the 


l 2 
minimum length criteria and see that these arbitrary components must be identically zero 


to attain the best approximation. 


— — Q0 0 0 
oO; on by 
b 1 
That is, x*e= f+} =]/0 — 0 0}} (9) 
C, C, bs 
0 0 0 0 0 
ba 
0 0 0 O 0 


The minimum length solution to Ax = pis, x” , Strang [3].Again, the useful result is 
specifically the equation that reveals x” . This process displays the matrix, which yields 


the desired result, 


I b 
— 0 00 ke 
co 0 0 0 6, o, 
0 a, 0 0 l b 
IfA= : then A* =| 0 — 0 Oland x*=A*b=|— (10) 
0 0 0 om C, 
0 000 0 0 00 0 
0 0 00 0 
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A’ is referred to as the pseudoinverse and is the matrix which provides for solution to the 


nearly singular system of Ax =b , Strang, [3]. 


This entire process has one sticking point. L* is the pseudoinverse described 


earlier. Now, the magnitude of rounding error in applying the inverse process is limited to 
the sum of the errors when inverting the individual siaoouel elements, o,. Each o,is the 
square root of each nonzero eigenvalue /, from both AA* and ATA . Now, the fine 
point, what happens when o, is sufficiently small to induce a divide by zero condition 
when taking the pseudoinverse? The reciprocal of o, is set to to zero by the code. Press, 
et. al., [7], denote this pecsdate as editing the singular values. The logic is sound. 

Recall the original formulation of the linear system. When redundant solutions (singular 
conditions) are encountered as a result of variable correlation, the matrix is unable to 
distinguish between the different basis functions and the associated distribution of the 
input data. By setting the reciprocal of any sufficiently close to zero singular value to 
zero, we effectively add.a zero multiple to the fitting parameters as opposed to some large 
combination of the basis functions that are degenerate to the best fit, [7]. Further, if any 


“nonzero singular value is very small, its reciprocal should also be set to zero. This term is 


most likely residual from rounding error and detracts from the optimal solution. 


A tule for determining the editing tolerance of singular values is given by Press, 


et. al., [7] as follows. Set to zero any singular value, o, when, 
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oO. 
The ratio LL. < Ne 


max 





where WN is the column length dimension and ¢ is the machine precision. 
Consider the following example least squares problem which illustrates the 
mathematical principles of the algorithm. In this example A is assumed symmetric on 


input just as in the Differential Correction form. 


3 -2 2 1 
Lett A=|-2 4 0O| and b=}{1 
2 QO 2 1 
where 
17 -14 10 
AfA=|-14 20 -4 
10 -4 8 


The eigenvalues 


A = 36=6° =o =6 
At A-AI are A, = 953 oo, =3 
A =0=0°=0, =() 


3 


whose corresponding eigenvectors are 


ih = ] 2 
v=-|-2|/,v =--|2!|,v =-} 1 

1 3 2 3 3 3 
] 2 —2 


Now A= U2V’or AV=U2 so 
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u =—Av, u =—Av 
1 ey H 2 Oo 2 
j 2 
3 -2 2|| 2 12 2 
u=-|-2 4 O}-2 eee 12|=-—|-2 
18 3 
2 O 2] 1 6 ] 
Similarly 
3 -2 211 3 ; l 
u =--|-2 4 OF; 2 ee =——|2 
2 9 3 
2 O 24/2 6 2 


Now for the zero singular value we must solve the homogeneous problem 


A’u =0 for u,. 
3 


l -2 24 
3 -2 2) 0 3 : 0 | 
—-2 4 0;0;>5)/0 1 Fi 0 | > 4, is arbitrary 
Z Qe Zey70 0 0 o| 2 
but, u= TU and u=-u 
—] 2 
1 } 
SO, u =u“,| --—|=—| 1 
3 3 5 3 
] ae 
Now lets refit the components. 
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T 
6 0 OT", 
A=UzVv’ =| u a |]o 3 Olw’ = bu 3usi«éOd = | 6u v! 3u v’ ov" 
l 2 3 2 1] 1 1 2 2 3 
0 0 Oll_r 


eae ee Mi ile. | 


This is simply the sum of rank one matrices. 


Now, 

2 ci 
3 6 2 -2 
a 2016 Ole “ms &. | 2 

A=Urv'= |—* —* 3 3 3 11 
Peale “ 
1 -2 2. 3 3 
a 8 


Now for the least squares piece. 
Ax=b and A=UEV' so. EV'x=U"b or x= VE~!U'b 
= | is the pseuodoinverse. Its structure was detailed earlier. 


Now applying equation (8) to the original right hand side yields, 


2 6 = 2 
9 9/1 9 

Oo “2 3 
Pia [7 
9 9 6 18 


This 1s the least squares solution of minimal length. 
Now the test, if x is orthogonal to M( A’), i.e., really the minimal length solution 


and u, is an element of M(A’ ) then orthogonality demands that their dot product equal 


Zero. 
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2 

3 

] 8 6 14 
— |= —+ 
31 54 54 54 
=e 
3 


xeu, = @ 


od 
co | WWI ol 


This example is easily extended to matrices of greater dimension. 
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IV. COMPARISON OF METHODS 


A. GENERAL 


The procedures used for computational comparison were NAVSPACECOMs 
algorithms AMA03, AMA04, and AMA06 coupled with a process driver program and the 
standard SVD algorithms as taken from Numerical Recipes by Press, et. al., [7]. During 
research of the actual NAVSPACECOM code, it was found that Numerical Recipes was 
referenced repeatedly throughout. This observation set the precedence for applying the 
basic SVD codes of Numerical Recipes for computational analysis. This common ground 


should serve to standardize any resulting code development. 


It should be noted that throughout chapter 15, Numerical Recipes [7] the SVD is 
recommended for least squares problems. Citing the following, “...solution of a least 
squares problem directly from the normal equations 1s rather susceptible to roundoff 
error.” And, In some applications, the normal equations are perfectly acceptable for 
linear least squares problems. However, in many cases the normal equations are very close 
to singular.” The authors detail with significant correlation the precise difficulties that 
generally arise in Differential Correction. They go on further detailing the exact 
complication the routine AMAO06 attempts to eliminate. The following excerpt from 


Numerical Recipes sums up the exact nature of this entire analysis: 
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A zero pivot element may be encountered during the solution of the normal 
equations, with Gauss-Jordan, in which you get no solution at all. Or, a very small pivot 


may occur in which case you typically get fitted parameters @, with very large 


magnitudes that are delicately balanced to cancel out almost precisely when the fitted 
function is evaluated. 


Why does this commonly occur? The reason is that, data do not clearly 
distinguish between two or more basis functions provided. If two such functions or two 
different combinations of functions happen to fit the data equally well - or badly — then 
the matrix A, is unable to distinguish between them and becomes singular. There is 
irony in the fact that least squares problems are both overdetermined (number of data 
points greater than the number of parameters) and underdetermined (ambiguous 
combinations of parameters exist). The ambiguities can be extremely hard to notice a 
priori in complicated problems. 


The SVD gives exactly what we need. In the overdetermined system SVD 
produces a solution that is the best approximation in the least squares sense. In the case 


of the underdetermined system, SVD produces a solution whose values (the @, ’s) are 


the smallest in the least squares sense also what we want. When some combination of 
the basis functions is irrelevant to the fit, that combination is driven down to a small, 
innocuous value, rather than pushed up to delicately canceling infinities. 


B. SENSITIVITY ANALYSIS 


The following sections present information that details the inner workings of the 
individual numerical solvers. The concept to focus on deals with the aspect of matrix 


sensitivity. The following simplified example details the problem explicitly (see Nyhoff 


[16]). 
; 2 6 x 8 
Ax =b is = 
; cai! a 


] 
The solution is obviously x= a . Now perturb the nght hand side very slightly. 
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— 2 6 x 8 
Ax = b 1S = 
fs sansa * Legian 


, a 10 ee 
Now the not so obvious solutionis x= | where the very small perturbation in the 


right hand side has altered the solution on the order of 10’ times the constant term of the 


second equation. 


In performing rounding error analysis, we note that computers operate in floating 
point arithmetic. That is, they convert every problem into a “nearby” problem perform 
calculations and then convert back the answer. The symbol ¢, called machine epsilon 
refers to the computer’s ability to distinguish between two consecutive binary 
representations of actual input values. When the relative representation of two 
consecutive numbers exceeds the computer’s ability to distinguish between them, 


( || a-d]| < €) we say an underflow has occurred in the calculation sequence. 


The following method, known as matrix perturbation theory, helps to analyze the 
nature of error in order to minimize it through algorithm refinement. Errors are 
accumulated primarily from two sources. First there may be measurement error sontaines 
within the original input data, and second, the algorithm itself may cause error through its 
internal approximations while processing calculations. For excellent derivations of the 


most commonly required error measurement techniques, see Demmel [6]. 
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The measure of the accumulated error is referred to as the condition number of a 
matrix. Simply put, the condition number indicates the precise level of sensitivity a matrix 


has relative to these types of very small perturbations. 


To summarize, an example taken from Strang [3] is in order. It directly illustrates 
the short fall in the present NAVSPACECOM code methodology of scaling the input 


matrix by the value 25,000,000. 


Begin by consider the linear system Ax = b , now perturb the nght hand side by 
Ob. These errors might have come from the observational data or roundoff. The change 


is small but the direction of change can not be controlled. The solution is subsequently 


changed from x to x + 6x. Now our system has become Aédx = db. We now must 
estimate the resulting perturbation 6x = Ab. There will be a large change in the 


solution when A”‘is large, i.e. A is nearly singular. Now consider our symmetric matrix 


with positive eigenvalues where0</A, <A, <---< A,. Any vector 6bis a combination of 


the corresponding unit eigenvectors x,,---,x, . Let ¢ indicate a very small change. 
ob 
If ob = &X, then wan % 


The error in |b] is amplified by rt which is the largest eigenvalue of A7'. The 


amplification is greatest when J, 1s close to zero. Thus nearly singular matrices are the 
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most sensitive. The serious drawback with this measure of sensitvity appears when scaling 
is introduced. Multiplying the matrix by a large scalar also scales the eigenvalue by the 
corresponding amount. This makes the matrix appear much less singular. This rescaling 


can’t however make an ill-conditioned matrix well. It is true dx will be smaller by the 


(ex 
x! 


same. The factor ||x||in the denominator normalizes the problem against such trivial 


scale factor however so will the solution to x = A™'b. The relative error “—~ stays the 


rescaling. There is a corresponding normalization of ||db||. The problem is to compare the 


relative change ——~* a | with the relative error ——* [6 a 


[bj Ix] 


At this point, I refer to the following theorem as taken from Strang, [3]. 
Theorem 


For a positive definite matrix, the solution x = A™'b and the error 5x = A™'dbalways 


satisfy, 


j= PE ants 
Ay, A, 


Therefore the relative error is bounded by 
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bx] _ 4, [6b 
Ik] 7 [bl 
Ae. Hise 2 sad 
The ratio c= ra = a is called the condition number of A. 


This analysis can be applied to the example given previously. 
C. ERROR ANALYSIS 


The following analysis applies to the least squares system. The established error 
bounds hold regardless of solution method. A fundamental rounding error analysis 
follows directly from the previous sensitivity analysis. See Golub, [5] and Demmel, [6] for 


very complete presentations of both rounding and backward error analysis 


The tradition definition of a vector or matrix norm holds throughout the following 
analysis. See Strang [3] as required. Where the symbol |||| appears in equation form, 


consistent use of the chosen norm is implied throughout. Condition number 1s as defined 


previously. 


Given the linear system Ax = b where r=b—Ax . Where A is nxn 


nonsingular, x is the computed solution, and 7 is the residual. Letx, be the exact solution 


and e the error. 


Then r=b-Ax= Ax, — Ax= A(x, ~ X) 
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or r= Ae 


Given A nonsingular then e = Ar. Taking consistent norms of the equations yields the 


following inequalities. 
lel<[A7Ip] and [Alllel > Fl 


Relating these inequalities bounds the error as follows, 


A“ = fel = 


[Al 
Now the exact solutions form, 
x,= Ab and Ax,=b 


Where 












































x, < [Ap and |Allx,|| = [bl 
We can now cout the exact solution. 
A“libl| = > Ib 
opm = hale Bt 


Now divide inequality (10) through by and form 











x, 
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(12) 





JA“Irl Ae o 


xX 


é 





























Ww | and in the right denomination by 


|A7 ilb| the expression for relative error now becomes, 





a a _ 


tly May aI ~ fala 


Where lel is the relative error. 


€ 
































Now recalling that the condition number of the matrix A can also be expressed as 


| cond(A) = |A“| [Al , which is required for analysis of methods that do not return 


eignevalue information, equation (11) becomes, 


rl), 2 (br 
conan) Et x,|| * cond(A) [bl 


Ib) 
This formulation is convenient in that it allows for testing of the maximum and 




















minimum associated errors for any given output as initiated by the condition of the matrix 


system. For more specific information on relative error and on absolute error given 
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machine dependent performance information see Engeln-Mullges and Uhlig [7]. This 


analysis can be applied to the given test matrices. 
D. COMPUTATIONAL EFFICIENCY 


Flops, or floating point operations, are the measure of algorithm efficiency. While 
they do not necessarily indicate the actual processing time, as different computers do 
internal processing differently, they do give an excellent measure of the magnitude of total 


number of computer calculation required by an algorithm as related to the functional input. 


In general, flop counts are obtained by summing the number of arithmetic 
operations from the most deeply nested algorithm statements. As an example of the 
accounting notice, a dot product operation of length n involves n multiplications and n 


additions. It therefore requires 2n flops. For matrix multiplications, the general form is 
CG, j) = AGk)* BE, 7) + CGS) 
This process requires 2mnp flops where CE R"”, AER”, Be R””. 


Applying the above flop counting procedure, the following flop count estimates 


are given for each algorithm as related to the matrix dimension [5]. 
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ALGORITHM FLOP COUNT 


Gauss-Jordan 


SVD 
(U,2,V) 

SVD 
(z,V only) 


This does not take into account specialized storage techniques and factorizations 












Am’n +8mn? + 9n° 







4mn’? + 8n° 
Minimum requirement for least squares 





when dealing with specific subprocess forms that may be optimized. In this respect, these 
estimates would be iain worst case. In some instances the total flop count of a 
particular subform may be cut by nearly half the indicated estimate. While ffisse methods 
yield respectable approximations of the computational cost, they do not provide actual run 
time analysis. Specific algorithms can be made extremely efficient when optimized for a 


particular computer. 
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V. CONCLUSIONS 


A. SUMMARY OF FINDINGS 


Testing of the actual NAVSPACECOM code and the SVD code was 
accomplished through PC based Digital Visual Fortran 6.0. All NAVSPACECOM source 
code was compiled and linked using the visual Fortran development environment. As 
noted earlier SVD source codes were adopted directly from Numerical Recipes [6]. 
Driver programs for both software suites were ieelones as adoptations of the LAPACK 
[9] and Numerical Recipes source code for advanced linear algebra. Test matrices were | 
generated using MATLAB and placed in standard text format for file upload by the driver 
programs. All Fortran source code was version 77. The smatiaiabie system used was a 


Micron, 300 MHz, Pentium PI. 


The test matrices indicated at Appendices A and B were used to determine how 
well the two routines compared when attempting to solve highly singular and known 
singular input matrices. As expected, the NAVSPACECOM code returns the SING = 0 
error and does not continue further. The SVD algorithm on the other hand returns. 
Further, upon analysis of the factored structure, the reconstructed matrices are within one 


order of magnitude of the original input. The differences are due to roundoff error. 


While the results presented in the output of the SVD algorithm may not look so 


appealing at first glance, it is important to note,that the present Gauss-Jordan, 
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NAVSPACECOM routines return no output at all. This is much more encouraging. As a 
minimum, even with the ridiculously singular nearly impossible input matrices, the SVD 
routine successfully returned the best approximation as the desired solution. At this point 
the differential correction subroutine now possesses an excellent starting point to continue 
processing this observational input set. The nonlinear method may then still converge. 
Had the NAVSPACECOM routine been allowed to process this set, the entire set would 
have been discarded based solely on the coincidental fact that the input matrix derivation is 
singular beyond the computer’s ability to distinguish otherwise. The macro effect of this 
NAVSPACECOM code shortfall is that potentially accurately derived observational data 
is now discarded. This will certainly effect the long-term distribution of the accumulated 


error in the orbital element set. 
B. RECOMMENDATIONS 


The potential for improvement in NAVSPACECOMs differential correction 
process through integration of an SVD algorithm which replace existing subroutines 
AMA03, AMA04, and AMAO6 exists. The only way to verify the magnitude of 
improvement, however, is to implement a test routine. The evidence presented herein 


indicates that further study should be pursued. 


C. CONCLUSION 


The results of this study are not all inclusive. Only the worst possible input matrix 
conditions were modeled. The frequency of occurrence of this type of input would 
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certainly have to be taken into consideration when deciding whether to upgrade the 
present process. Differential Correction involves highly complex series of numerical 
routines each with its own influence on the over all solution process. My findings indicate 
that some improvement may occur through reduced processing time and the numerical 
error of the actual values returned in the update state vectors. Additionally, it is worth 
pointing out that many modern statistical regression-fitting packages have begun to use 


more sophisticated linear algebra solvers for similar reasons as those addressed herein. 


While at best the SVD algorithm is approximately 24 times more computationally 
expensive in flops, it has the advantage of not requiring cyclic preconditioning to start the 
solution process. The present NAVSPACECOM routine has the potential to cycle for 
significant seeds prior to achieving SING = 0 stop criteria. If significant numbers of 
orbital element sets cause the AMA06 routine to process several iis through before 
either failing or calling the solver, the relative difference in run time between the two 
methods may be very small. This aspect of the study will require actually run time 


information to determine the flop count delta. 


SVD solutions exhibit less error in each iterate as a result of their “best 
approximation”. It is possible that over time, the resulting thorough evaluation of the 
observational data by an SVD synthesizing algorithm could effectively decrease the error 
in actual known position information as measured by the difference in the orbital element 


set and the corresponding test data. 
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0.9688 


0.3557 


0.0490 


0.7553 


0.8948 


0.2861 


— 0.2512 


0.9327 


3.4538 


2.2680 


1.7824 





APPENDIX A 


Test Case 1. Singular A Matrix 


Constructed random uniform (0,1) 


0.1310 


0.9408 


0.7019 


0.8477 


0.2093 


0.4551 


0.0811 


0.8511 


0.5620 


0.3193 


0.3749 


0.8678 


0.3722 


0.0737 


0.1998 


0.0495 


Matrix A: 


0.7654 0.5979 


0.3375 
0.2120 
0.8116 
0.6335 
0.1799 
0.2255 


0.4911 


0.9492 
0.2888 
0.8888 
0.1016 


0.0653 


0.2343 


0.9331 


0.0631 0.7666 
0.2642 0.6661 
0.9995 | 0.1309 
0.2120 0.0954 
0.4984 0.0149 
0.2905 0.2882 
0.6728 0.8167 


0.9580 0.9855 


Symmetric, nxn matrix P=A‘A. 


2.2680 1.7824 2.6181 2.6412 1.9559 2.2782 


3.0953 


1.5425 


1.9052 2.7916 2.2445 1.9391 


1.5425 1.4978 1.6401 1.6543 1.0673 1.0142 
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Highly singular A matrix: col(4) = (col(1)+col(3))/2 and col(8) = (col(2)+col(7))/2 


0.4488 
0.8035 
0.4164 
0.4715 
0.1121 
0.3716 
0.4489 


0.9183 


2.2731 
2.5172 


1.2783 





2.6181 1.9052 1.6401 2.1291 2.1478 1.5116 1.7757 
2.6412 2.7916 1.6543 2.1478 3.0720 1.8867 2.3445 2.5680 
1.9559 2.2445 1.0673 1.5116 1.8867 2.8209 1.9602 2.1023 
2.2782 1.9391 1.0142 1.6462 2.3445 1.9602 2.7791 2.3591 
2.2731 2.5172 1.2783 1.7757 2.5680 2.1023 2.3591 2.4382 
Decomposition Matrices: 
Matrix U 
-0.405972 0.466532 -0.366483 0.180020 -0.532409 -0.026839 -0.006785 0.408157 
-0.387027 -0.208141 0.614596 mune nee -0.301571 -0.411911 0.408982 0.006625 
0.238695 0.389374 0.258549 0.159762 0.718607 -0.126578 -0.006750 0.408139 
-0.322333 0.427960 -0.053982 0.169871 0.093150 -0.075934 0.013911 -0.816422 
-0.404391 0.039450 0.240801 -0.441264 -0.010121 0.762960 -0.001406 0.000463 
-0.326910 -0.511282 -0.089709 0.724015 0.092756 0.301448 -0.000586 0.000171 
-0.346071 -0.284110 -0.593342 -0.385717 0.303011 -0.199784 0.408545 0.006797 
-0.366547 -0.246130 0.010611 -0.210450 0.000580 -0.308105 -0.815805 -0.013969 
Diagonal of Matrix ¥ 
16.942978 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
Matrix V-Transpose 
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-0.405972 
0.467637 
-0.370840 
0.169917 
-0.575737 
0.261994 
-0.223596 


0.000238 


Check product against original matrix: 


-0.387027 


-0.210013 


0.614798 


-0.013430 


-0.269208 


-0.244425 


-0.360038 


-0.407559 


Original P Matrix: 


3.453800 
2.268000 
1.782400 
2.618100 
2.641200 
1.955900 
2.278200 


2.273100 


Product Matrix = U*2*(V-Transpose): 


2.268000 


3.095300 


1.542500 


1.905200 


2.791600 


2.244500 


1.939100 


2.517200 


-0.238695 
0.388574 
0.254242 
0.165458 
0.669202 
0.394398 
-0.310865 


0.000561 


1.782400 
1.542500 
1.497800 
1.640100 
1.684300 
1.067300 
1.014200 


1.278300 


-0.322333 
0.428113 
-0.058314 
0.167672 
0.201179 
0.727229 
0.334620 


-0.000351 


2.618100 
1.905200 
1.640100 
2.129100 
2.147800 
1.511600 
1.646200 


1.775700 


-0.404391 
0.038740 
0.255496 
-0.432943 
-0.065735 
0.369304 
0.664486 


-0.001333 


2.641200 
2.791600 
1.654300 
2.147800 
3.072000 
1.886700 
2.344500 


2.568000 
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-0.326910 
-0.511046 
-0.115354 

0.720162 

0.073310 

0.160414 

0.262515 


-0.000496 


1.955900 
2.244500 
1.067300 
1.511600 
1.886700 
2.820900 
1.960200 


2.102300 





-0.346070 -0.366547 


-0.282277 


-0.581027 


-0.406542 


0.312403 


-0.054058 


-0.175497 


-0.407884 


2.278200 


1.939100 


1.014200 


1.646200 


2.344500 


1.960200 


2.779100 


2.359100 


-0.246149 
0.016870 
-0.210019 
0.021403 
-0.148875 
-0.265547 


0.817022 


2.273100 
2.517200 
1.278300 
1.775700 
2.568000 
2.102300 
2.359100 


2.438200 


2.792422 
2.662113 
1.641832 
2.217128 
2.781547 
2.248604 
2.380401 


2.521245 


2.662113 


2.537884 


1.565216 


2.113665 


2.651745 


2.143672 


2.269319 


2.403590 


1.641832 
1.565216 
0.965332 
1.303582 
1.635438 
1.322089 
1.399580 


1.482391 


2.217128 
2.113665 
1.303582 
1.760356 
2.208493 
1.785347 
1.889991 


2.001818 


2.781547 


2.651745 


1.635438 


2.208493 
2.770714 
2.239847 
2.371130 


2.511426 


44 


2.248604 
2.143672 
1.322089 
1.785347 
2.239847 
1.810694 
1.916823 


2.030238 


2.380400 
2.269318 
1.399580 
1.88999] 
2.371130 
1.916823 
2.029173 


2.149235 





2.521245 . 


2.403590 


1.482391] 


2.001818 


2.511426 


2.030238 


2.149235 


2.276402 








1.000000 


0.900000 


0.090000 


0.009000 


0.000900 


0.000090 


0.000009 


0.000001 


Singular P matrix: col(1) = col(2)+col(3) + ... + col(8) 


APPENDIX B 


Test Case 2. Extreme Singularity 


0.900000 0.090000 0.009000 0.000900 


0.900000 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 
0.090000 
0.000000 
enoee 
0.000000 
0.000000 


0.000000 


0.000000 
0.000000 
0.009000 
0.000000 
0.000000 
0.000000 


0.000000 


0.000000 
0.000000 
0.000000 
0.000900 
0.000000 


0.000000 


0.000000 


0.000090 


0.000000 


0.000000 


0.000000 


0.000000 


0.000090 


0.000000 


0.000000 


0.000009 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 


0.000009 


0.000000 





0.000001 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 
0.000000 


0.000001 


Decomposition Matrices: 


Matrix U 


-0.726830 -0.298545 0.377237 0.221507 -0.175290 0.150237 -0.127672 0.348738 





-0.685806 
~-0.037087 


-0.003546 


0.302705 


0.334128 


-0.443863 


0.811978 


-0.223180 


-0.190945 


0.175290 


0.175288 


0.174926 


-0.150237 


-0.150237 


-0.150235 


0.127672 
0.127672 


0.127672 


-0.348738 
-0.348738 


-0.348738 


-0.836417 0.031208 -0.329137 


-0.000353 -0.089017 -0.019792 0.864991 0.288233 -0.149928 0.127670 -0.348738 


45 


-0.000035 -0.008940 


-0.000004 -0.000894 


-0.002062 


-0.000207 


0.000000 -0,000099 -0.000023 


Diagonal of Matrix 2 


1.860190 0.000000 


Matrix V-Transpose 
-0.676744 -0.731343 
0.634569 -0.521247 
0.187067 -0.224736 
0.322910 -0.377922 
0.008663 -0.010371 
0.000000 0.000017 
0.000000 -0.000001 


0.000000 0.000000 


Check product against original matrix: 


Original Matrix P: 
1.000000 0.900000 
0.900000 0.900000 
0.090000 0.000000 


0.009000 0.000000 


0.000000 


-0.084492 
-0.570616 
0.406687 
0.708234 
0.016396 
0.000174 
-0.000015 


0.000001] 


0.090000 
0.000000 
0.090000 


0.000000 


0.089682 
0.008994 


0.001000 


0.000000 


-0.003917 
-0.004328 
0.865145 
-0.499791 
-0.041196 
0.001730 
-0.000147 


0.000014 


0.009000 
0.000000 
0.000000 


0.009000 


-0.886388 


-0.091411 


-0.010180 


0.000000 


-0.000497 


-0.001723 


0.024989 


-0.038919 


0.998458 


0.030609 


-0.001470 


0.000138 


0.000900 


0.000000 


0.000000 


0.000000 


46 


-0.261431 


0.899086 


0.102642 


0.000000 


-0.000049 


~-0.000168 


0.002324 


-0.001936 


0.030457 


-0.999175 


-0.026647 


0.001378 


0.000090 


0.000000 


0.000000 


0.000000 


0.127375 
0.248332 


0.916846 


0.000000 


-0.000005 
-0.000017 
0.000231 
-0.000172 
0.002268 
-0.026540 
0.999247 


0.028215 


0.000009 
0.000000 
0.000000 


0.000000 


-0.348738 
-0.348622 


-0.385685 


0.000000 


-0.000001 
-0.000002 
0.000026 
-0.000019 
0.000243 
-0.002122 
0.028168 


-0.999601 


0.000001 
0.000000 
0.000000 


0.000000 








0.000900 


0.000090 


0.000009 


0.000001 


Product U*2*(V-Transpose): 


0.914987 


0.863342 


0.046687 


0.004464 


0.000444 


0.000044 


0.000004 


0.000000 


0.000000 
0.000000 
0.000000 


0.000000 


0.988808 
0.932996 
0.050454 
0.004824 
0.000480 
0.000048 
0.000005 


0.000001 


0.000000 


0.000000 


0.000000 


0.000000 


0.114237 


0.107789 


0.005829 


0.000557 


0.000055 


0.000006 


0.000001 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 


0.005295 
0.004997 
0.000270 
0.000026 
0.000003 
0.000000 
0.000000 


0.000000 


0.000900 


0.000000 


0.000000 


0.000000 


0.000672 


0.000634 


0.000034 


0.000003 


0.000000 


0.000000 


0.000000 


0.000000 


47 


0.000000 


0.000090 


0.000000 


0.000000 


0.000067 


0.000063 


0.000003 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 


0.000000 
0.000000 
0.000009 


0.000000 


0.000007 
0.000006 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 


0.000000 





0.000000 
0.000000 
0.000000 


0.000001 


0.000001 
0.000001 
0.000000 
0.000000 
0.000000 | 
0.000000 
0.000000 


0.000000 
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